Propagation of Uncertainty in Rigid Body Attitude Flows 



Taeyoung Lee*^, Nalin A. Chaturvedi^, Amit Sanyal, Melvin Leok*, and N. Harris McClamroch^ 



o 
o 

<N 
(N 



a 



> 

CO 
(N 

o 
o 

a 



X 



Abstract — Motivated by attitude control and attitude esti- 
mation problems for a rigid body, computational methods are 
proposed to propagate uncertainties in the angular velocity and 
the attitude. The nonlinear attitude flow is determined by Euler- 
Poincare equations that describe the rotational dynamics of the 
rigid body acting under the influence of an attitude dependent 
potential and by a reconstruction equation that describes the 
kinematics expressed in terms of an orthogonal matrix repre- 
senting the rigid body attitude. Uncertainties in the angular 
velocity and attitude are described in terms of ellipsoidal sets 
that are propagated through this highly nonlinear attitude flow. 
Computational methods are proposed, one method based on a 
local linearization of the attitude flow and two methods based 
on propagation of a small (unscented) sample selected from 
the initial uncertainty ellipsoid. Each of these computational 
methods is constructed using the Lie group variational integra- 
tor algorithm, viewed as a discretization of the attitude flow 
dynamics. Computational results are obtained that indicate (1) 
the strongly nonlinear attitude flow characteristics and (2) the 
limitations of each of these methods, and indeed any method, 
in providing effective global bounds on the nonlinear attitude 
flow. 

I. Introduction 

As an integrable system, the attitude dynamics of the free 
rigid body are reasonably well understood. However, if there 
is an attitude dependent potential that influences the rigid 
body, then the dynamics can be very complex. In this paper, 
such attitude dynamics are studied for a rigid body with an 
inertially fixed pivot acting under the influence of uniform 
and constant gravity; this model is subsequently referred to 
as the 3D pendulum [1]. 

A rigid 3D pendulum is a rigid body supported by a 
fixed, frictionless pivot, acted on by gravitational forces. The 
supporting pivot allows three degrees of rotational freedom 
of the pendulum. Uniform, constant gravity is assumed. 
The terminology 3D pendulum refers to the fact that the 
pendulum is a rigid body with three spatial dimensions and 
the pendulum has three rotational degrees of freedom. 

Two reference frames are introduced. An inertial reference 
frame has its origin at the pivot; the first two axes lie in the 
horizontal plane and the third axis is vertical in the direction 
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of gravity. A reference frame fixed to the pendulum body is 
also introduced. The origin of this body-fixed frame is also 
located at the pivot. In the body fixed frame, the moment of 
inertia of the pendulum, J G R^^^, is constant. This moment 
of inertia can be computed from the traditional moment 
of inertia of a translated frame whose origin is located at 
the center of mass of the pendulum using the parallel axis 
theorem. Since the origin of the body fixed frame is located 
at the pivot, principal axes with respect to this frame can be 
defined for which the moment of inertia is a diagonal matrix. 
Note that the center of mass of the 3D pendulum may or may 
not lie on one of the principal axes defined in this way. 

The dynamics of the 3D pendulum are given by the Euler 
equation that includes the moment due to gravity: 



JQ = JVt X 1] + mgp X R^es. 
The rotational kinematics equations are 
R = RS{n). 



(1) 



(2) 



Equations ([B and Q define the full dynamics of a rigid 
pendulum on the tangent bundle TS0{3). The attitude of the 
pendulum is represented by a rotation matrix R G SO (3), 
and the angular velocity in the body fixed frame is denoted 
by ne M^. 

In the above equations, the mass of the pendulum and 
the gravitational acceleration are denoted by m and g, 
respectively, and p G is the vector from the pivot to 
the mass center of the pendulum expressed in the body fixed 
frame. The vector 63 = [0,0, 1]^ is the direction of gravity 
in the inertial frame, so that R^es is the direction of gravity 
in the pendulum-fixed frame. For a given vector a G M^, the 
skew- symmetric matrix S{a) is defined as 
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so that S{a)b = a x 6 for a 6 G M^. 

There are two disjoint equilibria when the direction of 
gravity in the body fixed frame is collinear with the vector 
p. We define 



H = \ {R,n) eTS0{3)\R^es 
I = i{R,Q) eTS0{3)\R^es- 



P 



\\p\\ 



as hanging equilibria and inverted equilibria, respectively. 

The objective of this paper is to study the nonlinear 
attitude flow of the 3D pendulum dynamics by characterizing 



the propagation of uncertainty in the attitude and angu- 
lar velocity. Simple analytical bounds can be obtained for 
small, localized uncertainty over short time periods, but such 
bounds provide no insight into the global attitude dynamics. 
Uncertainty propagation can also be studied, within a prob- 
abilistic framework, using the Liouville partial differential 
equation, a special case of the Fokker-Planck equation, but 
this leads to significant solution difficulties. Our approach 
in this paper is to make use of computational tools, based 
on the recently introduced Lie group variational integrator 
for the 3D pendulum and on a framework for propagating 
suitably defined uncertainty ellipsoids [2]. We show that 
these computational tools are useful in some cases, but 
they have important limitations that arise from the complex 
attitude dynamics of the 3D pendulum. 

This line of research was motivated by our prior work 
on attitude control and estimation [3], [4]. The results in 
this paper demonstrate the complex global dynamics that can 
occur for an uncontrolled 3D pendulum. In this way, the 
challenges of attitude control and estimation are made clear, 
at least in the case that global results are desired. 

II. Problem formulation 
An uncertainty ellipsoid on TS0{3) is defined by 

S{R, n,P) = { {R, Q) e TS0{3) I x^p-^x < l} , 

where x = [(; SQ] e and ( = \og{FFR), SQ = Q - Q e 
R^. The center of the ellipsoid is given by the rotation matrix 
R and the angular velocity Cl; the matrix P G MP^^ is the 
uncertainty matrix that characterizes the size of the ellipsoid. 
In particular, if the initial attitude and angular velocity are 
known to lie in an initial uncertainty ellipsoid, we seek 
computational methods to propagate the uncertainty so that 
we obtain another uncertainty ellipsoid at a later time in 
which the resulting attitude and angular velocity of the 3D 
pendulum are expected to lie. 

III. Uncertainty propagation 

We study three methods to propagate the uncertainty 
set over the time interval [0,T]. All three methods make 
use of the Lie group variational integrator to propagate 
the angular velocity and the attitude [2]. The Lie group 
variational integrator is described by the following discrete 
update equations, 

hSiJQk + ^Mk) = FkJd - JdF^, (4) 

Rk-\-l = RkFki (5) 

m^+i = F^JQk + ^F^Mk + ^Mfc+i, (6) 

where the subscript k denotes the kth discrete variable for 
a fixed integration step size h e and F^ G 5*0 (3) 
is the relative attitude between two adjacent integration 
steps. The nonstandard moment inertia matrix is given by 
Jd = |tr[J]/3x3 - J G R^^^, and the moment due to 
the gravity is denoted by Mk = mgp x R^es G M^. For 
a given {Rk.^k), & is solved to find Fk G S0{3). Then 



(^fc+i, ^fc+i) is obtained by (0) and Q. This yields a map 
{Rkj^k) ^ {Rk-\-ij^k-\-i) and this process is repeated. 

A. Uncertainty propagation using linearization on [0, T] 

In this method, the uncertainty ellipsoid is propagated 
by updating both the center of the ellipsoid and the the 
uncertainty matrix P. The center of the uncertainty ellipsoid, 
denoted by i^^, is propagated according to the Lie group 
variational integrator, initialized by the center of the initial 
uncertainty ellipsoid. We then assume that the uncertainty el- 
lipsoid is sufficiently small that the flow of points within the 
uncertainty ellipsoid is well- approximated by the linearized 
flow of the Lie group variational integrator about this center 
solution, which we denote by 

where Xk = [Ck]SQk] ^ R^ and the matrix Ak G R^^^ 
depends on Rk^^k- The expression for A^ can be found 
in [5]. Using this linearization, the uncertainty matrix is 
propagated according to 

Pfc+i = AkPkAl (7) 

with initial condition given by the initial uncertainty matrix. 
In this way, the uncertainty ellipsoid is propagated and 
determined at time T. 

B. Uncertainty propagation using unscented method on 
[0,T] 

In this method, we compute the propagated uncertainty 
ellipsoid by using the Lie group variational integrator to 
propagate a small sample of points selected from the initial 
uncertainty ellipsoid. Following the conventional wisdom, we 
choose the 12 points corresponding to the intersection of the 
boundary of the initial uncertainty ellipsoid and its principal 
axes. This choice is informed by the fact that if the initial 
uncertainty ellipsoid is propagated by a linear flow, it will 
remain an ellipsoid, and it will coincide with the minimal 
volume ellipsoid containing the the 12 intersection points 
propagated by the same linear flow. 

More explicitly, suppose that the initial uncertainty ellip- 
soid is given by 

f(i?o,^o,Po). 

Let X' e R and (j)' = [(j)^^; (j)}^] G R^ be the i-th eigen- 
value and eigenvector, respectively, of the uncertainty matrix 
Po e R^""^ for i G {1,2..., 6}. The intersection of the 
corresponding principal axis and the ellipsoid boundary is 
Then, the 12 intersection points are given by, 

{ (i?o exp(±v^4), no ± Vx^^h) } i G {1, 2, . . . , 6} . 

We propagate these 12 initial conditions using the Lie group 
variational integrator to determine the 12 values of the 
attitude and angular velocity at time T. We then determine 
a minimal volume ellipsoid that contains all 12 values of 
the attitude and angular velocity at time T [6]. From this 
computed ellipsoid, the center attitude, the center angular 
velocity, and the uncertainty matrix at time T can be com- 
puted. 



C. Uncertainty propagation using unscented method with re- 
sampling on [0, T] 

The unscented method propagates the same sampled points 
obtained from the initial uncertainty elHpsoid throughout the 
entire time period [0,T]. In this modification, we partition 
the time period [0, T] into subintervals and we use the Lie 
group variational integrator to propagate the 12 sample points 
throughout each subinterval. At the end of each subinterval, 
we determine a minimal volume ellipsoid that contains all 12 
values of the attitude and angular velocity at the end of that 
subinterval. We then select a new set of sample points located 
on the principal axes of this new uncertainty ellipsoid. On 
the first subinterval, the 12 initial values of the attitude and 
angular velocity are obtained from the initial uncertainty 
ellipsoid. This method differs from the previous method in 
that resampling of the propagated 12 solutions is carried out 
at the end of each subinterval. 

It is worth remarking that if the flow is linear, the re- 
sampling technique will yield the same estimates as the usual 
unscented method. However, if the flow is nonlinear, and 
the re- sampling is sufficiently frequent, then the re- sampling 
technique will tend to yield a more conservative estimate. 
Indeed, the difference between the two estimates gives an 
indication of how nonlinear the flow is. 

IV. Numerical examples 

We apply these methods to the attitude dynamics of the 
3D pendulum. The pendulum body is chosen as an elliptic 
cylinder, and its properties are given by 



J = diag[0.13, 0.28, 0.17] kg • 



1 kg, p = O.Ses m. 



We consider two initial conditions; the first initial condition 
results in a near-oscillatory attitude flow while the second 
initial condition results in a highly irregular attitude flow. 

Oscillatory attitude flow: The initial uncertainty el- 
lipsoid is defined by its initial center attitude and angular 
velocity and the initial uncertainty matrix, which are given 
by 



^0 = ^3x3 

Po diag 



Qo = [3.0,0.1,0.1] rad/s. 



(5 — )^[1,1,1],0.01^[1,1,1] 

It is convenient for graphical purposes to consider the re- 
duced attitude R^es on the two-sphere and to plot the 
induced attitude flow on S^. This is due to the fact that the 3D 
pendulum has a symmetry of §^ given by a rotation about the 
vertical axis, and the configuration space can be reduced to a 
quotient space 5'0(3)/S^ using this symmetry [7]. The 
reduced attitude denotes the direction of gravity in the body 
fixed frame. Plots of the reduced attitude and the angular 
velocity responses, corresponding to the initial conditions, 
are shown in Fig. [T] for a time period of 10 seconds. 

We apply the three methods to propagate the initial uncer- 
tainty through this oscillatory attitude flow for 10 seconds. 
The integration step size is h = 0.005, corresponding to 
2000 time steps over the length of the simulation. For the 
unscented method, we find the minimal volume covering 




(a) Reduced attitude F 



(b) Angular velocity Q 



Fig. 1. Reduced attitude and angular velocity for oscillatory attitude flow 
(The center of the sphere is the hanging equilibrium when R^es = es.) 




(a) Magnitude of uncertainty matrix (b) Percentage of sample points con- 
tained in the computed uncertainty (%) 

Fig. 2. Uncertainty propagation through oscillatory attitude flow (Lin- 
earization: dotted, Unscented: blue, Unscented with re-sampling: red) 



ellipsoid every 0.1 seconds. For the unscented method with 
re-sampling, we find the minimal volume covering ellipsoid 
every 0.1 seconds, and we choose new set of sample points 
every 2.0 seconds. To provide a baseline for comparing the 
performance of each method, we choose sample points in 
the interior of the initial uncertainty ellipsoid. In particular, 
we choose 144 points on the level set of 
and we numerically integrate each of them. 

The magnitudes of the uncertainty matrix and the percent- 
age of sample points (out of 144) contained in the computed 
uncertainty ellipsoids are shown in Fig. O The propagation 
of uncertainties in the reduced attitude on are shown in 
Fig.E 

For a short time period, the properties of all methods are 
similar. The uncertainty ellipsoid computed by the unscented 
method is larger than for the linearization method, which is 
a reflection of the nonlinear nature of the dynamics. The 
re- sampling in the unscented method with re- sampling has 
the effect of enlarging the uncertainty ellipsoids. Thus, the 
magnitude of the uncertainty ellipsoid increases rapidly, and 
the uncertainty ellipsoid for the re-sampling method contains 
a greater proportion of the sample points after 6 seconds. 

The mean precentages of sample points contained in 
the computed uncertainty ellipsoids are 28.9%, 39.8%, and 
59.8%, respectively. This suggests that even for this oscilla- 
tory attitude flow, the nonlinear effects are so strong that it 
is difficult to accurately propagate the uncertainty. 

The computation times are 3.24, 45.18, and 45.84 seconds, 
respectively. 
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(a) Reduced attitude F 



(b) Angular velocity Q 



Fig. 3. Reduced attitude and angular velocity for irregular attitude flow 
(The center of the sphere is the hanging equilibrium when R^es = es.) 




(a) Magnitude of uncertainty matrix (b) Percentage of sample points con- 
tained in the computed uncertainty (%) 

Fig. 4. Uncertainty propagation through irregular attitude flow (Lineariza- 
tion: dotted, Unscented: blue, Unscented with re-sampling: red) 



Irregular attitude flow: The initial uncertainty ellipsoid 
is defined by its initial attitude and angular velocity and the 
initial uncertainty matrix, which are given by 



Po- 



diag (5 



TT 

180 



= [4.14,4.14,4.14] rad/s, 
)^[l,l,l],0.0l2[l,l,l] 



The reduced attitude on and the angular velocity responses 
for 10 seconds corresponding to the center initial conditions 
are shown in Fig. [S] 

We apply the three methods to propagate the initial un- 
certainty through this irregular attitude flow for 5 seconds. 
The integration step size is h = 0.002, corresponding to 
2500 time steps over the length of the simulation. For the 
unscented method, we find the minimal volume covering 
ellipsoid every 0.1 seconds. For the unscented method with 
re-sampling, we find the minimal volume covering ellipsoid 
every 0.1 seconds, and we choose new set of sample points 
every 0.5 seconds. To compare the properties of each method, 
we choose 144 initial sample points on the level set of 
= 0.8, and we numerically integrate them. 

The magnitudes of the uncertainty matrix and the percent- 
age of sample points (out of 144) contained in the computed 
uncertainty ellipsoid are shown in Fig. (H The propagation 
of uncertainties in the reduced attitude on are shown in 
Fig.H 

For a short time period, the properties of all methods 
are similar. The uncertainty ellipsoid computed by the lin- 
earization method grows rapidly, but it encloses few points 



after 1.5 second. The unscented method encloses more points 
with smaller uncertainty ellipsoids than for the linearization 
method, but it encloses few points after 3.5 second. The 
re- sampling in the unscented method with re- sampling has 
the effect of enlarging the uncertainty ellipsoids. Thus, the 
size of the uncertainty ellipsoids increase rapidly, and the 
uncertainty ellipsoid contains more sample points than the 
other two methods. 

The mean numbers of sample points contained in the com- 
puted uncertainty ellipsoid are 10.7%, 26.3%, and 73.55%, 
respectively. The computation times are 4.97, 43.50, and 
45.56 seconds, respectively. 

V. Global Features of the Attitude Flow 

The numerical results presented in the previous sec- 
tion demonstrate the difficulty in obtaining accurate global 
bounds on attitude solutions that are initialized in an un- 
certainty ellipsoid. It is claimed that the source of this 
difficulty is the nonlinear attitude flow of the 3D pendulum, 
especially the fact that the flow can exhibit chaos and extreme 
sensitivity to initial conditions. A conceptual description of 
certain global features of the attitude flow that help to explain 
this effect is now provided. 

As described in [7], the global dynamics of the 3D pen- 
dulum are complicated. There is a ID hanging equilibrium 
submanifold of the 3D configuration manifold, consisting 
of hanging equilibria that differ by a rotation about the 
vertical. There is also a ID inverted equilibrium submanifold 
consisting of inverted equilibria that differ by a rotation 
about the vertical. Each hanging equilibrium is stable in the 
sense of Liapunov. Each inverted equilibrium is unstable, 
with a 2D stable manifold, a 2D unstable manifold, and a 
2D center manifold. Let M denote the union of all the 2D 
stable manifolds corresponding to inverted equlibria. This 3D 
set M plays an important role in understanding the global 
dynamics of the attitude flow. 

Every trajectory in M converges to the inverted equi- 
librium manifold. Although the set M has zero measure, 
its existence influences the dynamics of the 3D pendulum 
attitude flow near M. Since M is constructed as the union 
of the stable manifolds of unstable equilibria, trajectories 
near M remain near M for an extended period of time. In 
particular, the closer a trajectory is to M the longer it remains 
near M. In fact, there are trajectories that remain close to 
M for arbitrarily long time periods. This property is due to 
the saddle character of each inverted equilibrium. 

It should be mentioned that it is difficult to determine 
exactly the set M. One can make use of linear attitude 
equations near an inverted equilibrium to approximate the 
tangent space to the stable manifold of that equilibrium. 
However, this provides only local information about M; 
the non-local properties of the set M are not understood. 
In practice, to accurately compute the global structure of a 
stable manifold, one relies on either (i) extremely high-order 
Taylor approximations of the nonlinear stable manifold for 
a neighborhood of the equilibrium, which is used to obtain 
sample points on the stable manifold that are then propagated 



backwards in time in order to compute the global structure of 
the stable manifold [8], or (ii) set-oriented techniques based 
on representing the nonlinear flow map for short times as a 
Markov process [9]. 

This argument demonstrates that the set M has a strong 
influence on the 3D pendulum dynamics near M, with 
high shearing and thus high sensitivity of the attitude flow 
near M. This is one of the mechanisms leading to the 
complex nonlinear dynamics of the 3D pendulum and makes 
it impossible to efficiently compute accurate global bounds 
on attitude solutions that are initialized in an uncertainty 
ellipsoid. 



VI. Conclusions 

The Lie group variational integrator is known to provide 
accurate long-term solutions of the rigid body equations in 
the presence of an external potential for a given initial atti- 
tude and angular velocity; these computed solutions exactly 
conserve the theoretical conservation properties, namely the 
symplectic structure, and the angular momentum component 
about the vertical axis in the case for the 3D pendulum. In 
addition, it exhibits very good energy behavior, with only 
a very small bounded energy oscillation, for exponentially 
long times. Furthermore, the Lie group variational integrator 
is also known to exactly conserve orthogonality of the 
computed attitude as a rotation matrix. 

Consequently, inaccuracies in the approximate ellipsoidal 
bounds computed according to the three computational meth- 
ods introduced do not arise from computational difficulties 
with the Lie group variational integrator. Rather, the inac- 
curacies in the approximate ellipsoidal bounds arise from 
the fact that the attitude flow dynamics are highly nonlinear, 
with regions wherein the dynamics cannot be adequately 
approximated by linear dynamics. 

With these qualifications, it is clear that the unscented 
method with resampling is the most accurate of the three pro- 
posed methods in propagating the uncertainty. This method 
can provide a basis for analysis of control and estimation 
problems for attitude systems such as the 3D pendulum. 
For example, digital control operates open loop between 
sample times; the analysis in this paper demonstrates the 
importance of the choice of inter- sampling time in obtaining 
accurate propagation of the attitude flow between sample 
times. As another example, attitude estimation operates 
open loop between measurement times; the analysis in this 
paper demonstrates the importance of the choice of inter- 
measurement time in obtaining accurate propagation of the 
attitude flow between measurement times. 

The bottom line demonstrated by the development in this 
paper is that guaranteed global bounds for attitude dynamics 
defined by the 3D pendulum, or indeed for any attitude 
dynamics with a nontrivial potential, are not achievable. That 
is, there is no universal approach to global and uniform 
approximation of the attitude flow dynamics. 
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Unscented with re-sampling 



Fig. 5. Uncertainty projected onto the reduced attitude on for oscillatory attitude flow (The center of the sphere is the hanging equilibrium when 
R'^es = 63.) 




t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 
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Fig. 6. Uncertainty projected onto the reduced attitude on for irregular attitude flow (The center of the sphere is the hanging equilibrium when 
R^es = 63.) 



